library(CausalGAM)
library(Matching)

M <- 1000
n.vec <- c(250, 500, 1000)

output.mat <- matrix(NA, 1000*3*3*2*7, 8)

counter <- 1
for (n in n.vec){
  for (alpha.ind in 1:3){
    for (beta.ind in 1:2){
      cat("n =         ", n, "\n")
      cat("alpha.ind = ", alpha.ind, "\n")
      cat("beta.ind =  ", beta.ind, "\n\n")
      for (i in 1:M){
        filename <- paste("../../data/MCdata.n", n, ".alpha", alpha.ind,
                          ".beta", beta.ind,
                          ".rep", i, ".Rdata", sep="")
        
        load(filename)

        for (spec in 1:7){
          if (spec == 1){
            pscore.out <- gam(x~s(z1,z2), family=binomial(probit), data=mydata)
            attach(mydata, warn.conflicts=FALSE)
            out <- Match(Y=y, Tr=x, X=predict(pscore.out), estimand="ATE")
            detach(mydata)
          }
          else if (spec == 2){
            pscore.out <- gam(x~s(z1,z2)+s(z3),
                              family=binomial(probit), data=mydata)
            attach(mydata, warn.conflicts=FALSE)
            out <- Match(Y=y, Tr=x, X=predict(pscore.out), estimand="ATE")
            detach(mydata)
          }
          else if (spec == 3){
            pscore.out <- gam(x~s(z1), family=binomial(probit), data=mydata)
            attach(mydata, warn.conflicts=FALSE)
            out <- Match(Y=y, Tr=x, X=predict(pscore.out), estimand="ATE")
            detach(mydata)
          }
          else if (spec == 4){
            pscore.out <- gam(x~s(z2), family=binomial(probit), data=mydata)
            attach(mydata, warn.conflicts=FALSE)
            out <- Match(Y=y, Tr=x, X=predict(pscore.out), estimand="ATE")
            detach(mydata)
          }
          else if (spec == 5){
            pscore.out <- gam(x~s(z2), family=binomial(probit), data=mydata)
            attach(mydata, warn.conflicts=FALSE)
            out <- Match(Y=y, Tr=x, X=predict(pscore.out), estimand="ATE")
            detach(mydata)
          }
          else if (spec == 6){
            pscore.out <- gam(x~s(z2)+s(z3),
                              family=binomial(probit), data=mydata)
            attach(mydata, warn.conflicts=FALSE)
            out <- Match(Y=y, Tr=x, X=predict(pscore.out), estimand="ATE")
            detach(mydata)
          }
          else if (spec == 7){
            pscore.out <- gam(x~s(z2), family=binomial(probit), data=mydata)
            attach(mydata, warn.conflicts=FALSE)
            out <- Match(Y=y, Tr=x, X=predict(pscore.out), estimand="ATE")
            detach(mydata)
          }
          
          output.mat[counter,1] <- out$est
          output.mat[counter,2] <- out$se
          output.mat[counter,3] <- (abs(out$est - 5.00) / out$se) < qnorm(.975)
          output.mat[counter,4] <- n
          output.mat[counter,5] <- alpha.ind
          output.mat[counter,6] <- i
          output.mat[counter,7] <- spec
          output.mat[counter,8] <- beta.ind
          
          counter <- counter + 1
        } ## end spec loop


        
      }
    }
  }
}
colnames(output.mat) <- c("ATE.hat", "ATE.se", "cover", "n", "alpha.ind", "rep",
                          "spec", "beta.ind")
write.csv(output.mat, file="Matching.csv")
